<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dsyevr.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="DSYEVR.1"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a>( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
     $                   ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
     $                   IWORK, LIWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          JOBZ, RANGE, UPLO
      INTEGER            IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
      DOUBLE PRECISION   ABSTOL, VL, VU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            ISUPPZ( * ), IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DSYEVR.22"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a> computes selected eigenvalues and, optionally, eigenvectors
</span><span class="comment">*</span><span class="comment">  of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
</span><span class="comment">*</span><span class="comment">  selected by specifying either a range of values or a range of
</span><span class="comment">*</span><span class="comment">  indices for the desired eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DSYEVR.27"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a> first reduces the matrix A to tridiagonal form T with a call
</span><span class="comment">*</span><span class="comment">  to <a name="DSYTRD.28"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>.  Then, whenever possible, <a name="DSYEVR.28"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a> calls <a name="DSTEMR.28"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> to compute
</span><span class="comment">*</span><span class="comment">  the eigenspectrum using Relatively Robust Representations.  <a name="DSTEMR.29"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>
</span><span class="comment">*</span><span class="comment">  computes eigenvalues by the dqds algorithm, while orthogonal
</span><span class="comment">*</span><span class="comment">  eigenvectors are computed from various &quot;good&quot; L D L^T representations
</span><span class="comment">*</span><span class="comment">  (also known as Relatively Robust Representations). Gram-Schmidt
</span><span class="comment">*</span><span class="comment">  orthogonalization is avoided as far as possible. More specifically,
</span><span class="comment">*</span><span class="comment">  the various steps of the algorithm are as follows.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  For each unreduced block (submatrix) of T,
</span><span class="comment">*</span><span class="comment">     (a) Compute T - sigma I  = L D L^T, so that L and D
</span><span class="comment">*</span><span class="comment">         define all the wanted eigenvalues to high relative accuracy.
</span><span class="comment">*</span><span class="comment">         This means that small relative changes in the entries of D and L
</span><span class="comment">*</span><span class="comment">         cause only small relative changes in the eigenvalues and
</span><span class="comment">*</span><span class="comment">         eigenvectors. The standard (unfactored) representation of the
</span><span class="comment">*</span><span class="comment">         tridiagonal matrix T does not have this property in general.
</span><span class="comment">*</span><span class="comment">     (b) Compute the eigenvalues to suitable accuracy.
</span><span class="comment">*</span><span class="comment">         If the eigenvectors are desired, the algorithm attains full
</span><span class="comment">*</span><span class="comment">         accuracy of the computed eigenvalues only right before
</span><span class="comment">*</span><span class="comment">         the corresponding vectors have to be computed, see steps c) and d).
</span><span class="comment">*</span><span class="comment">     (c) For each cluster of close eigenvalues, select a new
</span><span class="comment">*</span><span class="comment">         shift close to the cluster, find a new factorization, and refine
</span><span class="comment">*</span><span class="comment">         the shifted eigenvalues to suitable accuracy.
</span><span class="comment">*</span><span class="comment">     (d) For each eigenvalue with a large enough relative separation compute
</span><span class="comment">*</span><span class="comment">         the corresponding eigenvector by forming a rank revealing twisted
</span><span class="comment">*</span><span class="comment">         factorization. Go back to (c) for any clusters that remain.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The desired accuracy of the output can be specified by the input
</span><span class="comment">*</span><span class="comment">  parameter ABSTOL.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  For more details, see <a name="DSTEMR.57"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>'s documentation and:
</span><span class="comment">*</span><span class="comment">  - Inderjit S. Dhillon and Beresford N. Parlett: &quot;Multiple representations
</span><span class="comment">*</span><span class="comment">    to compute orthogonal eigenvectors of symmetric tridiagonal matrices,&quot;
</span><span class="comment">*</span><span class="comment">    Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
</span><span class="comment">*</span><span class="comment">  - Inderjit Dhillon and Beresford Parlett: &quot;Orthogonal Eigenvectors and
</span><span class="comment">*</span><span class="comment">    Relative Gaps,&quot; SIAM Journal on Matrix Analysis and Applications, Vol. 25,
</span><span class="comment">*</span><span class="comment">    2004.  Also LAPACK Working Note 154.
</span><span class="comment">*</span><span class="comment">  - Inderjit Dhillon: &quot;A new O(n^2) algorithm for the symmetric
</span><span class="comment">*</span><span class="comment">    tridiagonal eigenvalue/eigenvector problem&quot;,
</span><span class="comment">*</span><span class="comment">    Computer Science Division Technical Report No. UCB/CSD-97-971,
</span><span class="comment">*</span><span class="comment">    UC Berkeley, May 1997.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Note 1 : <a name="DSYEVR.70"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a> calls <a name="DSTEMR.70"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> when the full spectrum is requested
</span><span class="comment">*</span><span class="comment">  on machines which conform to the ieee-754 floating point standard.
</span><span class="comment">*</span><span class="comment">  <a name="DSYEVR.72"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a> calls <a name="DSTEBZ.72"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and <a name="SSTEIN.72"></a><a href="sstein.f.html#SSTEIN.1">SSTEIN</a> on non-ieee machines and
</span><span class="comment">*</span><span class="comment">  when partial spectrum requests are made.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Normal execution of <a name="DSTEMR.75"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> may create NaNs and infinities and
</span><span class="comment">*</span><span class="comment">  hence may abort due to a floating point exception in environments
</span><span class="comment">*</span><span class="comment">  which do not handle NaNs and infinities in the ieee standard default
</span><span class="comment">*</span><span class="comment">  manner.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  JOBZ    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'N':  Compute eigenvalues only;
</span><span class="comment">*</span><span class="comment">          = 'V':  Compute eigenvalues and eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RANGE   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'A': all eigenvalues will be found.
</span><span class="comment">*</span><span class="comment">          = 'V': all eigenvalues in the half-open interval (VL,VU]
</span><span class="comment">*</span><span class="comment">                 will be found.
</span><span class="comment">*</span><span class="comment">          = 'I': the IL-th through IU-th eigenvalues will be found.
</span><span class="comment">*</span><span class="comment">********* For RANGE = 'V' or 'I' and IU - IL &lt; N - 1, <a name="DSTEBZ.92"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and
</span><span class="comment">*</span><span class="comment">********* <a name="DSTEIN.93"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> are called
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  UPLO    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'U':  Upper triangle of A is stored;
</span><span class="comment">*</span><span class="comment">          = 'L':  Lower triangle of A is stored.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix A.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
</span><span class="comment">*</span><span class="comment">          On entry, the symmetric matrix A.  If UPLO = 'U', the
</span><span class="comment">*</span><span class="comment">          leading N-by-N upper triangular part of A contains the
</span><span class="comment">*</span><span class="comment">          upper triangular part of the matrix A.  If UPLO = 'L',
</span><span class="comment">*</span><span class="comment">          the leading N-by-N lower triangular part of A contains
</span><span class="comment">*</span><span class="comment">          the lower triangular part of the matrix A.
</span><span class="comment">*</span><span class="comment">          On exit, the lower triangle (if UPLO='L') or the upper
</span><span class="comment">*</span><span class="comment">          triangle (if UPLO='U') of A, including the diagonal, is
</span><span class="comment">*</span><span class="comment">          destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A.  LDA &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VL      (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">  VU      (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          If RANGE='V', the lower and upper bounds of the interval to
</span><span class="comment">*</span><span class="comment">          be searched for eigenvalues. VL &lt; VU.
</span><span class="comment">*</span><span class="comment">          Not referenced if RANGE = 'A' or 'I'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IL      (input) INTEGER
</span><span class="comment">*</span><span class="comment">  IU      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          If RANGE='I', the indices (in ascending order) of the
</span><span class="comment">*</span><span class="comment">          smallest and largest eigenvalues to be returned.
</span><span class="comment">*</span><span class="comment">          1 &lt;= IL &lt;= IU &lt;= N, if N &gt; 0; IL = 1 and IU = 0 if N = 0.
</span><span class="comment">*</span><span class="comment">          Not referenced if RANGE = 'A' or 'V'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ABSTOL  (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The absolute error tolerance for the eigenvalues.
</span><span class="comment">*</span><span class="comment">          An approximate eigenvalue is accepted as converged
</span><span class="comment">*</span><span class="comment">          when it is determined to lie in an interval [a,b]
</span><span class="comment">*</span><span class="comment">          of width less than or equal to
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                  ABSTOL + EPS *   max( |a|,|b| ) ,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          where EPS is the machine precision.  If ABSTOL is less than
</span><span class="comment">*</span><span class="comment">          or equal to zero, then  EPS*|T|  will be used in its place,
</span><span class="comment">*</span><span class="comment">          where |T| is the 1-norm of the tridiagonal matrix obtained
</span><span class="comment">*</span><span class="comment">          by reducing A to tridiagonal form.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          See &quot;Computing Small Singular Values of Bidiagonal Matrices
</span><span class="comment">*</span><span class="comment">          with Guaranteed High Relative Accuracy,&quot; by Demmel and
</span><span class="comment">*</span><span class="comment">          Kahan, LAPACK Working Note #3.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If high relative accuracy is important, set ABSTOL to
</span><span class="comment">*</span><span class="comment">          <a name="DLAMCH.146"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( 'Safe minimum' ).  Doing so will guarantee that
</span><span class="comment">*</span><span class="comment">          eigenvalues are computed to high relative accuracy when
</span><span class="comment">*</span><span class="comment">          possible in future releases.  The current code does not
</span><span class="comment">*</span><span class="comment">          make any guarantees about high relative accuracy, but
</span><span class="comment">*</span><span class="comment">          future releases will. See J. Barlow and J. Demmel,
</span><span class="comment">*</span><span class="comment">          &quot;Computing Accurate Eigensystems of Scaled Diagonally
</span><span class="comment">*</span><span class="comment">          Dominant Matrices&quot;, LAPACK Working Note #7, for a discussion
</span><span class="comment">*</span><span class="comment">          of which matrices define their eigenvalues to high relative
</span><span class="comment">*</span><span class="comment">          accuracy.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The total number of eigenvalues found.  0 &lt;= M &lt;= N.
</span><span class="comment">*</span><span class="comment">          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The first M elements contain the selected eigenvalues in
</span><span class="comment">*</span><span class="comment">          ascending order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M))
</span><span class="comment">*</span><span class="comment">          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
</span><span class="comment">*</span><span class="comment">          contain the orthonormal eigenvectors of the matrix A
</span><span class="comment">*</span><span class="comment">          corresponding to the selected eigenvalues, with the i-th
</span><span class="comment">*</span><span class="comment">          column of Z holding the eigenvector associated with W(i).
</span><span class="comment">*</span><span class="comment">          If JOBZ = 'N', then Z is not referenced.
</span><span class="comment">*</span><span class="comment">          Note: the user must ensure that at least max(1,M) columns are
</span><span class="comment">*</span><span class="comment">          supplied in the array Z; if RANGE = 'V', the exact value of M
</span><span class="comment">*</span><span class="comment">          is not known in advance and an upper bound must be used.
</span><span class="comment">*</span><span class="comment">          Supplying N columns is always safe.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDZ     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array Z.  LDZ &gt;= 1, and if
</span><span class="comment">*</span><span class="comment">          JOBZ = 'V', LDZ &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) )
</span><span class="comment">*</span><span class="comment">          The support of the eigenvectors in Z, i.e., the indices
</span><span class="comment">*</span><span class="comment">          indicating the nonzero elements in Z. The i-th eigenvector
</span><span class="comment">*</span><span class="comment">          is nonzero only in elements ISUPPZ( 2*i-1 ) through
</span><span class="comment">*</span><span class="comment">          ISUPPZ( 2*i ).
</span><span class="comment">*</span><span class="comment">********* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK.  LWORK &gt;= max(1,26*N).
</span><span class="comment">*</span><span class="comment">          For optimal efficiency, LWORK &gt;= (NB+6)*N,
</span><span class="comment">*</span><span class="comment">          where NB is the max of the blocksize for <a name="DSYTRD.192"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a> and <a name="DORMTR.192"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>
</span><span class="comment">*</span><span class="comment">          returned by <a name="ILAENV.193"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.198"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, IWORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LIWORK  (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array IWORK.  LIWORK &gt;= max(1,10*N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LIWORK = -1, then a workspace query is assumed; the
</span><span class="comment">*</span><span class="comment">          routine only calculates the optimal size of the IWORK array,
</span><span class="comment">*</span><span class="comment">          returns this value as the first entry of the IWORK array, and
</span><span class="comment">*</span><span class="comment">          no error message related to LIWORK is issued by <a name="XERBLA.209"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  Internal error
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Inderjit Dhillon, IBM Almaden, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">     Ken Stanley, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment">       California at Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Jason Riedy, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment">       California at Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 0.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
     $                   TRYRAC
      CHARACTER          ORDER
      INTEGER            I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
     $                   INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
     $                   INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
     $                   LLWORK, LLWRKN, LWKOPT, LWMIN, NB, NSPLIT
      DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
     $                   SIGMA, SMLNUM, TMP1, VLL, VUU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.245"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.246"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.247"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANSY.247"></a><a href="dlansy.f.html#DLANSY.1">DLANSY</a>
      EXTERNAL           <a name="LSAME.248"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.248"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>, <a name="DLAMCH.248"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANSY.248"></a><a href="dlansy.f.html#DLANSY.1">DLANSY</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY, <a name="DORMTR.251"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>, DSCAL, <a name="DSTEBZ.251"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>, <a name="DSTEMR.251"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>, <a name="DSTEIN.251"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>,
     $                   <a name="DSTERF.252"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>, DSWAP, <a name="DSYTRD.252"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>, <a name="XERBLA.252"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, MIN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      IEEEOK = <a name="ILAENV.261"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 10, <span class="string">'<a name="DSYEVR.261"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a>'</span>, <span class="string">'N'</span>, 1, 2, 3, 4 )
<span class="comment">*</span><span class="comment">
</span>      LOWER = <a name="LSAME.263"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> )
      WANTZ = <a name="LSAME.264"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBZ, <span class="string">'V'</span> )
      ALLEIG = <a name="LSAME.265"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> )
      VALEIG = <a name="LSAME.266"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> )
      INDEIG = <a name="LSAME.267"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> )
<span class="comment">*</span><span class="comment">
</span>      LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
<span class="comment">*</span><span class="comment">
</span>      LWMIN = MAX( 1, 26*N )
      LIWMIN = MAX( 1, 10*N )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.( WANTZ .OR. <a name="LSAME.275"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBZ, <span class="string">'N'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( LOWER .OR. <a name="LSAME.279"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) ) ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE
         IF( VALEIG ) THEN
            IF( N.GT.0 .AND. VU.LE.VL )
     $         INFO = -8
         ELSE IF( INDEIG ) THEN
            IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
               INFO = -9
            ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
               INFO = -10
            END IF
         END IF
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
            INFO = -15
         ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -18
         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -20
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         NB = <a name="ILAENV.308"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DSYTRD.308"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>'</span>, UPLO, N, -1, -1, -1 )
         NB = MAX( NB, <a name="ILAENV.309"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DORMTR.309"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>'</span>, UPLO, N, -1, -1, -1 ) )
         LWKOPT = MAX( ( NB+1 )*N, LWMIN )
         WORK( 1 ) = LWKOPT
         IWORK( 1 ) = LIWMIN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.316"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DSYEVR.316"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      IF( N.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         WORK( 1 ) = 7
         IF( ALLEIG .OR. INDEIG ) THEN
            M = 1
            W( 1 ) = A( 1, 1 )
         ELSE
            IF( VL.LT.A( 1, 1 ) .AND. VU.GE.A( 1, 1 ) ) THEN
               M = 1
               W( 1 ) = A( 1, 1 )
            END IF
         END IF
         IF( WANTZ )
     $      Z( 1, 1 ) = ONE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants.
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="DLAMCH.348"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      EPS = <a name="DLAMCH.349"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SMLNUM = SAFMIN / EPS
      BIGNUM = ONE / SMLNUM
      RMIN = SQRT( SMLNUM )
      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale matrix to allowable range, if necessary.
</span><span class="comment">*</span><span class="comment">
</span>      ISCALE = 0
      ABSTLL = ABSTOL
      VLL = VL
      VUU = VU
      ANRM = <a name="DLANSY.361"></a><a href="dlansy.f.html#DLANSY.1">DLANSY</a>( <span class="string">'M'</span>, UPLO, N, A, LDA, WORK )
      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
         ISCALE = 1
         SIGMA = RMIN / ANRM
      ELSE IF( ANRM.GT.RMAX ) THEN
         ISCALE = 1
         SIGMA = RMAX / ANRM
      END IF
      IF( ISCALE.EQ.1 ) THEN
         IF( LOWER ) THEN
            DO 10 J = 1, N
               CALL DSCAL( N-J+1, SIGMA, A( J, J ), 1 )
   10       CONTINUE
         ELSE
            DO 20 J = 1, N
               CALL DSCAL( J, SIGMA, A( 1, J ), 1 )
   20       CONTINUE
         END IF
         IF( ABSTOL.GT.0 )
     $      ABSTLL = ABSTOL*SIGMA
         IF( VALEIG ) THEN
            VLL = VL*SIGMA
            VUU = VU*SIGMA
         END IF
      END IF

<span class="comment">*</span><span class="comment">     Initialize indices into workspaces.  Note: The IWORK indices are
</span><span class="comment">*</span><span class="comment">     used only if <a name="DSTERF.388"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="DSTEMR.388"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> fail.
</span>
<span class="comment">*</span><span class="comment">     WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
</span><span class="comment">*</span><span class="comment">     elementary reflectors used in <a name="DSYTRD.391"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>.
</span>      INDTAU = 1
<span class="comment">*</span><span class="comment">     WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
</span>      INDD = INDTAU + N
<span class="comment">*</span><span class="comment">     WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
</span><span class="comment">*</span><span class="comment">     tridiagonal matrix from <a name="DSYTRD.396"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>.
</span>      INDE = INDD + N
<span class="comment">*</span><span class="comment">     WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
</span><span class="comment">*</span><span class="comment">     -written by <a name="DSTEMR.399"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> (the <a name="DSTERF.399"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> path copies the diagonal to W).
</span>      INDDD = INDE + N
<span class="comment">*</span><span class="comment">     WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
</span><span class="comment">*</span><span class="comment">     -written while computing the eigenvalues in <a name="DSTERF.402"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> and <a name="DSTEMR.402"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>.
</span>      INDEE = INDDD + N
<span class="comment">*</span><span class="comment">     INDWK is the starting offset of the left-over workspace, and
</span><span class="comment">*</span><span class="comment">     LLWORK is the remaining workspace size.
</span>      INDWK = INDEE + N
      LLWORK = LWORK - INDWK + 1

<span class="comment">*</span><span class="comment">     IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in <a name="DSTEBZ.409"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and
</span><span class="comment">*</span><span class="comment">     stores the block indices of each of the M&lt;=N eigenvalues.
</span>      INDIBL = 1
<span class="comment">*</span><span class="comment">     IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in <a name="DSTEBZ.412"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and
</span><span class="comment">*</span><span class="comment">     stores the starting and finishing indices of each block.
</span>      INDISP = INDIBL + N
<span class="comment">*</span><span class="comment">     IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
</span><span class="comment">*</span><span class="comment">     that corresponding to eigenvectors that fail to converge in
</span><span class="comment">*</span><span class="comment">     <a name="DSTEIN.417"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.  This information is discarded; if any fail, the driver
</span><span class="comment">*</span><span class="comment">     returns INFO &gt; 0.
</span>      INDIFL = INDISP + N
<span class="comment">*</span><span class="comment">     INDIWO is the offset of the remaining integer workspace.
</span>      INDIWO = INDISP + N

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Call <a name="DSYTRD.424"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a> to reduce symmetric matrix to tridiagonal form.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DSYTRD.426"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>( UPLO, N, A, LDA, WORK( INDD ), WORK( INDE ),
     $             WORK( INDTAU ), WORK( INDWK ), LLWORK, IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If all eigenvalues are desired
</span><span class="comment">*</span><span class="comment">     then call <a name="DSTERF.430"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="DSTEMR.430"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> and <a name="DORMTR.430"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( ALLEIG .OR. ( INDEIG .AND. IL.EQ.1 .AND. IU.EQ.N ) ) .AND.
     $    IEEEOK.EQ.1 ) THEN
         IF( .NOT.WANTZ ) THEN
            CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
            CALL <a name="DSTERF.437"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>( N, W, WORK( INDEE ), INFO )
         ELSE
            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
            CALL DCOPY( N, WORK( INDD ), 1, WORK( INDDD ), 1 )
<span class="comment">*</span><span class="comment">
</span>            IF (ABSTOL .LE. TWO*N*EPS) THEN
               TRYRAC = .TRUE.
            ELSE
               TRYRAC = .FALSE.
            END IF
            CALL <a name="DSTEMR.447"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>( JOBZ, <span class="string">'A'</span>, N, WORK( INDDD ), WORK( INDEE ),
     $                   VL, VU, IL, IU, M, W, Z, LDZ, N, ISUPPZ,
     $                   TRYRAC, WORK( INDWK ), LWORK, IWORK, LIWORK,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply orthogonal matrix used in reduction to tridiagonal
</span><span class="comment">*</span><span class="comment">        form to eigenvectors returned by <a name="DSTEIN.455"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>            IF( WANTZ .AND. INFO.EQ.0 ) THEN
               INDWKN = INDE
               LLWRKN = LWORK - INDWKN + 1
               CALL <a name="DORMTR.460"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>( <span class="string">'L'</span>, UPLO, <span class="string">'N'</span>, N, M, A, LDA,
     $                      WORK( INDTAU ), Z, LDZ, WORK( INDWKN ),
     $                      LLWRKN, IINFO )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>         IF( INFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">           Everything worked.  Skip <a name="DSTEBZ.468"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>/<a name="DSTEIN.468"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.  IWORK(:) are
</span><span class="comment">*</span><span class="comment">           undefined.
</span>            M = N
            GO TO 30
         END IF
         INFO = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Otherwise, call <a name="DSTEBZ.476"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and, if eigenvectors are desired, <a name="DSTEIN.476"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">     Also call <a name="DSTEBZ.477"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and <a name="DSTEIN.477"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> if <a name="DSTEMR.477"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> fails.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         ORDER = <span class="string">'B'</span>
      ELSE
         ORDER = <span class="string">'E'</span>
      END IF

      CALL <a name="DSTEBZ.485"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
     $             WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
     $             IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWK ),
     $             IWORK( INDIWO ), INFO )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         CALL <a name="DSTEIN.491"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>( N, WORK( INDD ), WORK( INDE ), M, W,
     $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
     $                WORK( INDWK ), IWORK( INDIWO ), IWORK( INDIFL ),
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply orthogonal matrix used in reduction to tridiagonal
</span><span class="comment">*</span><span class="comment">        form to eigenvectors returned by <a name="DSTEIN.497"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>         INDWKN = INDE
         LLWRKN = LWORK - INDWKN + 1
         CALL <a name="DORMTR.501"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>( <span class="string">'L'</span>, UPLO, <span class="string">'N'</span>, N, M, A, LDA, WORK( INDTAU ), Z,
     $                LDZ, WORK( INDWKN ), LLWRKN, IINFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If matrix was scaled, then rescale eigenvalues appropriately.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Jump here if <a name="DSTEMR.507"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>/<a name="DSTEIN.507"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> succeeded.
</span>   30 CONTINUE
      IF( ISCALE.EQ.1 ) THEN
         IF( INFO.EQ.0 ) THEN
            IMAX = M
         ELSE
            IMAX = INFO - 1
         END IF
         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If eigenvalues are not in order, then sort them, along with
</span><span class="comment">*</span><span class="comment">     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
</span><span class="comment">*</span><span class="comment">     It may not be initialized (if <a name="DSTEMR.520"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>/<a name="DSTEIN.520"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> succeeded), and we do
</span><span class="comment">*</span><span class="comment">     not return this detailed information to the user.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         DO 50 J = 1, M - 1
            I = 0
            TMP1 = W( J )
            DO 40 JJ = J + 1, M
               IF( W( JJ ).LT.TMP1 ) THEN
                  I = JJ
                  TMP1 = W( JJ )
               END IF
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( I.NE.0 ) THEN
               W( I ) = W( J )
               W( J ) = TMP1
               CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
            END IF
   50    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set WORK(1) to optimal workspace size.
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWKOPT
      IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DSYEVR.549"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
